/*******************************************************
* REPLICATION CODE for main models in the article in Journal of European Public Policy:
* Mårtensson, M.; Palme, J.; Ruhs, M. and Österman, M. 'Shielding Free Movement? Reciprocity in Welfare Institutions and Opposition to EU Labour Immigration'
*******************************************************/
* This do-file reproduces the regression tables and the figures in the main article as well as the descriptive statistics presented in the online appendix.
* The tables are shown as Stata output and written to .rtf files and the figures as pdf and eps.
* Required packages
*	estout 
*	sencode
* NB: The code produce the tables in a raw format and the variable labels and the format do not correspond to the tables in the published article.

*NB: Since some temporary variables are used the complete code should be run at once

*set working directory where data is located (this is also where the figures and tables will be saved)
cd "?"

*Load dataset 
use "Replication_data_shielding_JEPP.dta", clear

*Control variables macro
global ind_controls "c.gndr c.agea c.agea#c.agea c.edlevdummy c.lrscale c.edctn_1 c.dsbld_1 c.rtrd_1 c.cmsrv_1 c.hswrk_1 c.dngoth_1 c.brncntr c.facntr c.mocntr c.self_emp"


***********************************************************************
*Descriptive statistics
***********************************************************************
*Country-level data descriptive statistics 

*Create temp VoC-variables for desc stat
tempvar MME
tempvar LME
tempvar CME
gen `MME'=(VoC==0)
gen `LME'=(VoC==1)
gen `CME'=(VoC==2)
collapse reciproindex_2010 ue_reciproindex_2010 gdpcapTEN_log  `LME' `CME' `MME', by(cntry)
eststo clear
estpost sum   reciproindex_2010 ue_reciproindex_2010 gdpcapTEN_log  `LME' `CME' `MME'
esttab using "descriptives_cntry.rtf", replace cells("count mean(fmt(a3))  sd(fmt(a3))  min(fmt(a3)) max(fmt(a3))")  noobs nomtitle nonumber 

*Indiviudal-level data descriptive statistics
use "Replication_data_shielding_JEPP.dta", clear
eststo clear
estpost sum frmvagg_2014 frmv_profdummy  uempl_joined_1 gndr agea edlevdummy lrscale pdwrk_1 edctn_1 dsbld_1 rtrd_1 cmsrv_1 hswrk_1 dngoth_1 self_emp brncntr facntr mocntr hincfel_01 if !mi(frmv_profdummy) & !mi(frmvagg_2014) & !mi(reciproindex_2010) [aweight=dweight]
esttab using "descriptives.rtf", replace cells("count mean(fmt(a3))  sd(fmt(a3))  min(fmt(a3)) max(fmt(a3))")  noobs nomtitle nonumber 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("frmvagg_2014") to("Opposition to free movement") 
sleep 300
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("reciproindex_2010") to("Reciprocity") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("ue_Reciproindex_2010") to("Reciprocity in unemployment insurance") 
sleep 300
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("frmv_profdummy") to("Professional dummy (experiment)") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("gdpcapTEN_log") to("GDP/cap logged") 
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("uempl_joined_1") to("Unemployed") 
sleep 300
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("gndr") to("Female") 
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("agea") to("Age") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("edlevdummy") to("Low education (dummy)") 
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("lrscale") to("Left-right scale") 
sleep 300
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("pdwrk_1") to("Main activity: paid work") 
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("edctn_1") to("Main activity: in education") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("dsbld_1") to("Main activity: permanently sick or disabled") 
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("rtrd_1") to("Main activity: retired") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("cmsrv_1") to("Main activity: community or military service") 
sleep 300
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("hswrk_1") to("Main activity: housework, looking after children") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("dngoth_1") to("Main activity: other") 
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("self_emp") to("Self employed") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("brncntr") to("Born in country") 
sleep 300
filefilter "descriptives1.rtf" "descriptives.rtf", replace from("facntr") to("Father born in country") 
filefilter "descriptives.rtf" "descriptives1.rtf", replace from("mocntr") to("Mother born in country") 
filefilter "descriptives1.rtf" "descriptives_individual.rtf", replace from("hincfel_01") to("Difficulty making ends meet on household income") 
erase descriptives.rtf 
erase descriptives1.rtf

***********************************************************************
*Main regression table, 
*Table 1. Regression analysis of opposition to EU labour immigration on institutional reciprocity in 12 EU and EFTA receiving countries
***********************************************************************
use "Replication_data_shielding_JEPP.dta", clear
eststo clear
eststo: reg frmvagg_2014 frmv_profdummy  reciproindex_2010   [pweight=dweight] , cluster(cntry)
estadd local cntry_FEs "No"
estadd local reci_IA "No"
estadd local ind_controls "No"

eststo: reg frmvagg_2014 frmv_profdummy  reciproindex_2010  gdpcapTEN_log  i.VoC [pweight=dweight] , cluster(cntry)
estadd local cntry_FEs "No"
estadd local reci_IA "No"
estadd local ind_controls "No"

eststo: reg frmvagg_2014 frmv_profdummy  reciproindex_2010 uempl_joined_1  gdpcapTEN_log i.VoC [pweight=dweight] , cluster(cntry)
estadd local cntry_FEs "No"
estadd local reci_IA "No"
estadd local ind_controls "No"

eststo: reg frmvagg_2014 frmv_profdummy  reciproindex_2010 uempl_joined_1  gdpcapTEN_log i.VoC  c.reciproindex_2010#c.uempl_joined_1 [pweight=dweight] , cluster(cntry)
estadd local cntry_FEs "No"
estadd local reci_IA "No"
estadd local ind_controls "No"

eststo: reg frmvagg_2014 frmv_profdummy reciproindex_2010 uempl_joined_1  gdpcapTEN_log i.VoC c.reciproindex_2010#c.uempl_joined_1 $ind_controls  [pweight=dweight] , cluster(cntry)
estadd local cntry_FEs "No"
estadd local reci_IA "No"
estadd local ind_controls "Yes"

eststo: reg frmvagg_2014 frmv_profdummy uempl_joined_1 c.reciproindex_2010#c.uempl_joined_1  $ind_controls i.cntry_num  [pweight=dweight] , cluster(cntry)
estadd local cntry_FEs "Yes"
estadd local reci_IA "No"
estadd local ind_controls "Yes"

eststo: reg frmvagg_2014 frmv_profdummy uempl_joined_1  c.reciproindex_2010#c.uempl_joined_1  $ind_controls c.reciproindex_2010#(c.frmv_profdummy $ind_controls)  i.cntry_num  [pweight=dweight] , cluster(cntry)
estadd local cntry_FEs "Yes"
estadd local reci_IA "Yes"
estadd local ind_controls "Yes"

*esttab word, for main text
esttab using "table_main.rtf", replace  se   stats(r2_a  N ind_controls cntry_FEs reci_IA, fmt(a2) labels("Adj. R$^2$" "Observations" "Individual-level controls" "Country FEs" "Reciprocity * ind. controls")) b(%9.2fc) ///
compress star(* 0.10 ** 0.05 *** 0.01) unstack label nomtitles  keep(_cons frmv_profdummy reciproindex_2010 uempl_joined_1 gdpcapTEN_log 1.VoC 2.VoC c.reciproindex_2010#c.uempl_joined_1) /// 
order(reciproindex_2010 gdpcapTEN_log 1.VoC 2.VoC frmv_profdummy uempl_joined_1  edlevdummy rtrd gndr agea lrscale c.reciproindex_2010#c.uempl_joined_1 ) ///
varlabels(c.reciproindex_2010#c.uempl_joined_1  "Reciprocity * UE")

*esttab word, for appendix with all coefficients
esttab using "table_main_app.rtf", replace  se   stats(r2_a  N ind_controls cntry_FEs reci_IA, fmt(a2) labels("Adj. R$^2$" "Observations" "Individual-level controls" "Country FEs" "Reciprocity * ind. controls")) b(%9.2fc) ///
compress star(* 0.10 ** 0.05 *** 0.01) unstack label nomtitles  keep(_cons frmv_profdummy reciproindex_2010 uempl_joined_1 gdpcapTEN_log 1.VoC 2.VoC c.reciproindex_2010#c.uempl_joined_1  ///
gndr agea c.agea#c.agea edlevdummy lrscale edctn_1 dsbld_1 rtrd_1 cmsrv_1 hswrk_1 dngoth_1 brncntr facntr mocntr self_emp) /// 
order(reciproindex_2010 gdpcapTEN_log 1.VoC 2.VoC frmv_profdummy uempl_joined_1  edlevdummy rtrd gndr agea lrscale c.reciproindex_2010#c.uempl_joined_1 ) ///
varlabels( c.reciproindex_2010#c.uempl_joined_1  "Reciprocity * UE")



*PLOTS
******************************************************************************
*Figure 1. Opposition to EU labour immigration (labourers and professionals) 
******************************************************************************
use "Replication_data_shielding_JEPP.dta", clear

*Labourers
*Create value for EU mean-level
qui su allbpe if essround==7 [aweight=dweight]
local EU_mean_lab = r(mean)

*Create numerical country identifier (for collapse)
encode cntry, gen (ncntry)

*Collapse mean, by country. 
collapse (mean) meanallbpe_2014=allbpe [pweight = dweight], by(ncntry)

*Add EU-mean as observation
set obs 13
replace ncntry=13 in 13
replace meanallbpe_2014 = `EU_mean_lab' in 13
label define ncntry 13 "Mean", add
label values ncntry ncntry

*Go back to string for country
decode ncntry, gen(cntry)
drop ncntry
 
*Sort
sort meanallbpe_2014
sencode cntry, replace

*Draw Plot
label define meanallbpe_2014 1 `" "1" "Many" "' 2 `" "2" "Some" "' 3 `" "3" " A few" "' 4 `" "4" "None "'
label values meanallbpe_2014 meanallbpe_2014
twoway (bar meanallbpe_2014 cntry if cntry !=7, bcolor(gs9) horizontal barwidth(0.70)) || ///
(bar meanallbpe_2014 cntry if cntry==7, bcolor(gs5) horizontal barwidth(0.70)), ytitle(Labourers, size(4)) ///
ylabel(#13, labsize(small) labels angle(horizontal) valuelabel noticks) ymtick(, noticks) xtitle(, size(zero)) ///
xlabel(1/4, labsize(small) valuelabel grid gmax) legend(off) scheme(s1mono) plotregion(margin(zero) fcolor(none) lcolor(none)) ///
saving(labourers, replace)
 
* Professionals
use "Replication_data_shielding_JEPP.dta", clear

*Create value for EU mean-level
qui su alpfpe if essround==7 [aweight=dweight]
local EU_mean_prof = r(mean)

*Create numerical country identifier (for collapse)
encode cntry, gen (ncntry)

*Collapse mean, by country. 
collapse (mean) meanalpfpe_2014=alpfpe [pweight = dweight], by(ncntry)

*Add EU-mean as observation
set obs 13
replace ncntry=13 in 13
replace meanalpfpe_2014 = `EU_mean_prof' in 13
label define ncntry 13 "Mean", add
label values ncntry ncntry

*Go back to string for country
decode ncntry, gen(cntry)
drop ncntry
 
*Sort
sort meanalpfpe_2014
sencode cntry, replace

*Draw Plot
label define meanalpfpe_2014 1 `" "1" "Many" "' 2 `" "2" "Some" "' 3 `" "3" " A few" "' 4 `" "4" "None "'
label values meanalpfpe_2014 meanalpfpe_2014
twoway (bar meanalpfpe_2014 cntry if cntry !=8, bcolor(gs9) horizontal barwidth(0.70)) || ///
(bar meanalpfpe_2014 cntry if cntry==8, bcolor(gs5) horizontal barwidth(0.70)), ytitle(Professionals, size(4)) ///
ylabel(#20, labsize(small) labels angle(horizontal) valuelabel noticks) ymtick(, noticks) xtitle(, size(zero)) ///
xlabel(1/4, labsize(small) valuelabel grid gmax) legend(off) scheme(s1mono) plotregion(margin(zero) fcolor(none) lcolor(none)) ///
saving(professionals, replace)
 
 *Combine plots
gr combine labourers.gph professionals.gph, cols(1) scheme(s1mono)
gr export figure1_opposition.pdf , replace
gr export figure1_opposition.eps, replace
 

*Scatterplot: opposition to free movement and reciprocity
****************************************************************
use "Replication_data_shielding_JEPP.dta", clear
set scheme plotplainblind
collapse frmvagg_2014 reciproindex_2010 if t==2, by(cntry)
scatter frmvagg_2014 reciproindex_2010 , mlabel(cntry) xlabel(0.3(0.1)0.8) ylabel(1(1)4) ytitle("Opposition to free movement") xtitle("General reciprocity") || lfit frmvagg_2014 reciproindex_2010, legend(off)
gr export "scatter_frmv_recip.pdf", replace
gr export "scatter_frmv_recip.eps", replace


*Marginal effects plot, general reciprocity, Model (7) in Table 1
******************************************************************
use "Replication_data_shielding_JEPP.dta", clear
*Create temporary variable for age squared
tempvar age_sq
gen `age_sq'= agea * agea

*Actual plot: run model + margins + marginsplot
global ind_controls "c.gndr c.agea c.agea#c.agea c.edlevdummy c.lrscale c.edctn_1 c.dsbld_1 c.rtrd_1 c.cmsrv_1 c.hswrk_1 c.dngoth_1 c.brncntr c.facntr c.mocntr c.self_emp"

reg frmvagg_2014 frmv_profdummy uempl_joined_1  c.reciproindex_2010#c.uempl_joined_1  $ind_controls c.reciproindex_2010#(c.frmv_profdummy $ind_controls)  i.cntry_num  [pweight=dweight] , cluster(cntry)

margins, dydx(uempl_joined_1) at(reciproindex_2010=(0.30(0.05)0.75))

marginsplot, yline(0) xlabel(0.3(0.1)0.8) title("") xtitle(Reciprocity) title(General reciprocity, size(medlarge)) addplot(hist reciproindex_2010, freq yaxis(2) yscale(alt axis(2) range(0 20000))  below  width(0.01) color(gs10) legend(off) ) name(main_ia, replace) 

*gr export "marginsplot, main recip.pdf", replace
*gr export "marginsplot, main recip.eps", replace


*Marginal effects plot, unemployment insurance reciprocity, Model (7) in Table A6 in online appendix
******************************************************************
 reg frmvagg_2014 frmv_profdummy uempl_joined_1  c.ue_reciproindex_2010#c.uempl_joined_1  $ind_controls c.ue_reciproindex_2010#(c.frmv_profdummy $ind_controls)  i.cntry_num  [pweight=dweight] , cluster(cntry)
margins, dydx(uempl_joined_1) at(ue_reciproindex_2010=(0.30(0.05)0.75))
marginsplot, yline(0) xlabel(0.3(0.1)0.8) title("") xtitle(Reciprocity UI) title(Reciprocity in unemployment insurance, size(medlarge)) addplot(hist ue_reciproindex_2010, freq yaxis(2) yscale(alt axis(2) range(0 20000))  below  width(0.01) color(gs10) legend(off) ) name(ue_ia, replace)  // ytitle("") ytitle(,axis(2) )
*gr export "marginsplot, ue recip.pdf", replace
*gr export "marginsplot, ue recip.eps", replace


*Combine marginal effect plots
**********************************************
*gr combine main_ia ue_ia,  xsize(7)
gr combine main_ia ue_ia, col(1) ysize(7)
gr export "marginsplot, gen ue comb.pdf", replace
gr export "marginsplot, gen ue comb.eps", replace


